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ABSTRACT 

Panel  methods  and  their  underlying  theory  are  reviewed  with 
regard  to  hydrodynamic  analysis  of  propeller  performance. 

Green's  identity  is  used  to  convert  the  differential  Laplace's 
equation  into  an  integral  equation.  The  velocity  potential  on 
the  surface  of  the  lifting  body  can  be  expressed  by 
integrating  the  potential  induced  by  source/doublet 
singularities  distributed  over  the  surface.  The  numerical 
discretizations  of  the  boundary  surface,  singularity 
distributions,  the  integral  equation,  and  the  formulation  of 
the  panel  method  are  discussed.  The  advantages  of  the 
application  of  panel  methods  in  viscous /inviscid  interactive 
procedure  and  propeller  blade  design  are  outlined.  Results  of 
propeller  blade  analysis  with  the  panel  method  are  presented, 
comparing  the  predictions  of  the  VSAERO  panel  method 
and  a  vortex  lattice  method  with  experimental  data.  The 
panel  method,  which  includes  consideration  of  propeller  hub 
effects,  gives  predictions  in  good  agreement  with 
experimental  data. 

ADMINISTRATIVE  INFORMATION 

This  investigation  was  sponsored  by  the  Chief  of  Naval  Research,  Office  of 
Naval  Technology  (Code  OCNR23)  under  the  Ship  and  Submarine  Technology 
Program,  Program  Element  62543N,  ONT  thrust  area  RS43-434  Propeller 
Quieting.  The  work  was  performed  at  the  David  Taylor  Research  Center  under 
Work  Unit  1508-001. 


INTRODUCTION 

Laplace’s  equation  is  one  of  the  most  frequently  encountered  equations  in 
the  field  of  engineering.  It  governs  the  potential  of  an  electrostatic  field,  the 
stress  function  of  torsion,  and  the  potential  of  an  incompressible  inviscid 
irrotational  flow  field.  The  methods  of  solution  have  been  known  for  quite 
some  time  and  can  be  formulated  in  either  differential  or  integral  forms.  The 
exact  solutions  can  be  found  for  a  number  of  problems  with  simplified 
geometries.  The  problem  of  interest  -  marine  propeller  hydrodynamics  -  involves 
extremely  complicated  geometries,  and  solutions  can  be  obtained  only 
numerically.  The  availability  of  modern  high-speed  computing  machinery  and 
the  maturity  of  some  numerical  techniques  have  resulted  in  the  recent 
development  of  panel  methods.1'5 

The  panel  method  has  a  distinct  advantage  over  the  differential  approaches, 
such  as  finite-element  or  finite-difference  methods,  because  the  unknowns  of  the 
panel  method  are  situated  only  on  the  fluid/solid  interface  and  not  throughout 
the  external  space.  In  principle,  the  method  applies  only  to  incompressible, 
inviscid,  and  irrotational  flow.  However,  due  to  recent  developments  of 
inviscid/viscid  interactive  techniques,  the  application  can  be  extended  to 
problems  with  mild  boundary  layer  separation  and  cavitation.  In  the  propeller 
hydrodynamics  application,  when  the  geometry  of  a  given  propeller  blade  and  its 
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advance  coefficient  are  specified,  the  panel  method  is  capable  of  calculating  the 
pressure,  lift,  drag,  and  moment  on  blades,6  duct,  or  band.  With  the 
implementation  of  a  wake  relaxation  procedure,  the  wake  contraction  can  also 
be  predicted.5  If  the  total  number  of  unknowns  can  be  kept  to  a  tolerable  level, 
the  panel  method  can  also  be  applied  to  the  problems  of  propeller/rudder  and 
propeller/hull  interactions. 

In  the  propeller  hydrodynamics  application,  as  reported6  earlier,  the  panel 
method  is  capable  of  providing  reliable  surface  pressure  distributions  in  the 
blade  leading-edge  region.  It  was  also  shown  that  the  method  is  capable  of 
predicting  the  leading -edge  pressure  peak  when  the  propeller  is  operating  at  off- 
design  advance  ratio,  although  comparison  with  experimental  data  is 
insufficient.  The  purpose  of  this  report  is  to  reinforce  Hess’s  observations 
through  comparison  of  various  solutions  obtained  by  the  panel  method,  the 
lifting  surface  theory,  the  equivalent  two-dimensional  theory,  and  extensive 
experimental  data  acquired  at  the  David  Taylor  Research  Center  (DTRC).  Some 
details  of  the  panel  method  formulation  are  outlined. 

MATHEMATICAL  BACKGROUND 
GOVERNING  EQUATIONS 

For  steady,  inviscid,  irrotational,  and  incompressible  flows  in  the  domain  Q 
bounded  by  the  surface  S,  there  exists  a  velocity  potential  that  satisfies  the 
Laplace’s  equation,  with  appropriate  boundary  conditions  at  S, 


V2  0  =  0 


(1) 


where 


V2_.1L  _1L  _JL 

dx2  +  9y2  +  3z2 

and  x,  y,  and  z  are  the  orthogonal  coordinate  axes.  The  assumptions  made  in 
deriving  Eq.  1  and  the  consequence  of  the  assumptions  are  described  in  the 
Appendix. 

BOUNDARY  CONDITIONS 

Two  types  of  boundary  conditions  are  often  found  in  hydrodynamic 
applications:  (a)  the  Neuman-type  boundary  condition  (specifying  the  derivative 
of  the  potential  O)  applied  at  the  Kutta  point  or  on  the  wetted  surface  and  (b) 
the  Dirichlet-type  boundary  condition  (specifying  the  potential)  applied  at  the 
internal  flow  to  render  a  unique  solution  of  external  flow  problems  via  Green’s 
third  identity.  The  Dirichlet-type  boundary  condition  may  also  arise  in  design 
procedures  in  which  the  surface  pressure  distribution  is  specified. 

The  idealized  flow  and  its  boundaries  can  be  represented  as  shown  in  Fig.  1. 
On  the  solid  surface,  the  velocity  component  normal  to  the  surface  is  set  to 
zero  for  a  nonpermeable  condition,  or  to  the  boundary  layer  transpiration 
velocity  if  inviscid/viscid  interaction  is  implemented.  Wakes  that  carry  away  the 
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vorticity  generated  in  the  boundary  layer  on  the  body  may  be  assumed  to  attach 
to  the  sharp  edges,  such  as  the  trailing  edge,  leading  edge,  tip,  or  a  line  of 
separation.  Suppose  in  Fig.  1  that  the  potential  is  Oj  in  the  external  region  <t>,, 

<t>2  in  the  internal  region  Q2,  and  <t>3  in  the  wake  region  Qr  The  potentials  <t>2 
and  <J>3  do  not  simulate  the  real  flow;  only  their  influences  on  the  external  region 
are  simulated. 


INTEGRAL  EQUATIONS 

Suppose  a  domain  Q  is  bounded  by  a  surface  S  as  shown  in  Fig.  2.  If  f  and 
g  are  any  two  similar  fields  which  possess  continuous  second  derivatives,  then 
Green’s  identity  states 

f  (gWf-fV2g)dQ  =  f  (ffrVg-gn-Vf)dS  (2) 

J<2  J  S 

with  the  unit  vector  If  normal  to  the  surface  S  and  pointing  into  the  domain  Q. 
Suppose  g  is  replaced  by  a  velocity  potential  0  that  satisfies  Laplace’s  equation 
V2<t>  =  0  and  f  is  replaced  by  1/r,  where  r  is  the  length  of  a  vector  "r  from  any 
source  point  Q  interior  to  the  domain  2  to  a  fixed  point  P.  The  function  1/r 
possesses  second-order  derivatives  and  satifies  Laplace’s  equation  V2(l / r)  =  0  at  all 
P,  except  when  P  is  interior  to  Q  or  on  S. 

When  P  is  exterior  to  Q,  straightforward  application  of  Eq.  2  results  in 

f  [1/r  n  V0-<t>n-V(l/r)]dS  =  0.  (3) 

Js 


When  P  is  interior  to  Q,  1/r  becomes  singular  as  Q  approaches  P.  The 
singularity  can  be  avoided  by  surrounding  the  fixed  point  P  with  a  small  sphere 
of  radius  t,  surface  S',  and  volume  Q',  and  applying  Eq.  2  to  the  volume 
Q-Q'  and  surface  S  +  S',  to  obtain 

f  [(1/r)  n  VO-<Dn-V(l/r)]dS  =  0.  (4) 

Js  +  S' 


As  the  radius  z  of  the  sphere  surrounding  P  approaches  zero,  the  following 
limits  can  be  obtained: 


and 


(1/c)  n-V$  dS' 


-  0 


(5a) 


<t> n-7  (1/r)  dS' 


-4n<t>(P). 


(5b) 


Equations  4  and  S  imply 


4rr<t>(P)  = 


[<Pn  V(l/r)  -  l/rn-V4>]dS. 


(6) 
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When  P  lies  on  the  surface  S,  1/r  becomes  singular  as  Q  approaches  P.  The 
singularity  can  be  avoided  in  a  manner  similar  to  that  just  described.  If  the 
surface  S  is  smooth  (first-order  derivatives  are  continuous),  the  following  limits 
can  be  obtained: 

f  (1/e)  n-VO  dS'  -  0  (7a) 

J  S' 

and 

f  <t>n-V(l/r)  dS'  -*•  ~2n<t>(P).  (7b) 

Js' 

It  is  concluded  that 

10:  P  lies  exterior  of  Q 

4n<t>(P):  P  lies  interior  of  Q 
2n<t>(P):  P  lies  on  S.  (8) 

The  function  f  in  Eq.  2  was  introduced  by  Green  in  solving  Laplace’s  equation 
V2g  =  0  and  was  called  Green’s  function  later  by  Riemann.  It  is  clear  from  Eq.  8 
that  with  the  help  of  Green’s  function  the  potential  at  any  fixed  point  P  can  be 
expressed  in  terms  of  its  value  on  the  surface.  Such  a  method  of  solving 
Laplace’s  equation,  as  opposed  to  the  method  of  using  series  of  special 
functions,  is  called  the  method  of  singularities. 

The  geometry  of  practical  problems  is  complicated  and  often  involves 
multiple  domains.  Figure  1  represents  the  geometry  of  a  lifting  body  problem.  It 
involves  three  distinct  domains,  Q,,  Q2  and  Qr  The  domain  Q,  is  the  external 
flow  field,  Qj  is  the  lifting  body  separated  from  52,  by  the  boundary  surface  S,2, 
and  Qj  is  the  wake  separated  from  Q,  by  the  surface  S13  and  from  Q2  by  the 
surface  S2J.  As  mentioned  previously,  the  potentials  <t>2  and  <t>3  in  Q2  and  Q3  are 
used  to  simulate  effects  on  the  external  flow  field  Q,  and  not  the  real  flow  of  a 
lifting  body  or  wake. 

Consider  a  fixed  point  P  lying  in  the  external  flow  field  S2,  and  apply  Eq.  8 
in  terms  of  the  potential  <t>,: 


4i*t>(P)  =  f  -  (l/r)  n, -70,  dS  + 

f  <D,n,-V(l/r)  dS 

^S12 

Js<2 

+  f  -  (l/r)  n,  V<t>,  dS  + 

f  <J>,n,  V(l/r)  dS 

JS,3 

^SI3 

+  f  -  (l/r)  n, -VO,  dS  + 
Jso 

f  4>,n,-V(l/r)  dS  (9) 

Jso 

Since  P  lies  in  Q,  and  is  outside  of  S2,  and  since, 
in  terms  of  the  potential  $2  gives 

n^=  -H,  on  S,2,  Eq.  8  applied 

4 


(1/r)  dS 


0=  -  f  -  (1/r)  n,  V«J)2  dS  +  f  -  ot,-V 

Js.2  Js12 

+  [  -  (1/r)  dS  +  f  V  (1/r)  dS 

Js.. 


(10) 


Since  P  is  outside  of  £23  and  n‘3=  -n^  on  S13,  if  we  apply  Eq.  8  in  terms  of  the 
potential  <t>3  we  have 


•V4>3  dS  + 


f  - 

J  S  . 


V  (1/r)  dS 
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o  =  -  f  -  (1/r)  n,' 

K 

+  J  -  (l/OiyVOj  dS  +  f  Ojlij-V  (1/r)  dS 

JS«  K 

Summing  the  contribution  of  all  surfaces  to  the  potential  at  point  P,  we  have 
4n<t>(P)  =  j  -  (1/r)  nf  ■  (VU),  -  V$2)  dS  +  f  («!>,  -  <t>2)  n,  •  V  (1/r)  dS 

K  K 

+  I  -  (1/r)  n,  •  (VO,  -  V03)  dS  +  f  (<t>,  -  <t>3)  n,  V  (1/r)  dS 

Js,3  JS,3 

+  f  -  (1/r)  n^-(V<t>2- V03)  dS  +  f  (<t>2  -  03)  iy  V  (1/r)  dS 

Js23  K 

+  f  -  (1/r)  n.-VO,  dS  +  f  <D  n,  V  (1/r)  dS. 

Jsn  Jsn 


(11) 


(12) 


If  the  outer  surface  S0  lies  infinitely  far  away  from  the  other  surfaces  and  the 
point  P  lies  near  Sfl,  the  contributions  of  surfaces  SI2,  SJ3,  and  S^  to  the  point  P 
become  negligible  and  Eq.  12  becomes 


(13) 


4i*t>w  =  f  -  (1/r)  nj-  70,  dS  +  f  d^n,*  V  (1/r)  dS 

J  S0~*«»  Jso~*-oo 

where  represents  the  unperturbed  potential  in  the  far  field.  If  the  wake  is 
assumed  to  be  infinitely  thin,  the  surface  length  S23  approaches  zero  and  its 
contribution  to  the  potential  at  point  P  vanishes.  The  wake  effect  on  the 
potential  at  point  P  is  represented  by  the  surface  integral  on  S13  in  Eq.  12: 

f  -  (l/r)tv(V<t>,-V<t>3)  dS  +  f  <•,-*,)  n  -v  (1/r)  dS.  (14) 

K  K 

Expression  14  can  be  decomposed  to  upper  (  +  )  and  lower(-)  surfaces: 
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f  (l/r)  n+  (V0(+  -V03+)  dS  +  f  (O;  -03+)  n+  -V  (1/r)  dS 

Jsr3  H 

+  f  (l/r)  n '  •  (VO- -  V03')  dS  +  f  (0;  - O")  n"  •  V  (l/r)  dS.  (15) 

Jsr3  jsr3 

In  the  limit  of  infinitesimal  wake  thickness,  the  following  situations  arise: 

n+  =  -la- 


<t>+  —  <h- 
3  v3  * 


Equation  15  then  becomes 


f  (l/r)n+(V0,+  -V10-)  dS+  +  f  (Of+  -0,-)n+  V  (l/r)  dS.  (16) 

Js13  Js13 

If  the  wakes  are  assumed  to  be  nonload-carrying  with  no  flow  entrainment, 
then  the  potential  gradient  across  the  wake  surface  is  continuous  and  the  first 
term  in  Eq.  16  vanishes.  The  potential  jump  across  the  wake  surface  in  the 
second  term  in  Eq.  16  can  be  replaced  by  AO  and  Eq.  12  can  then  be  written  as 

4nO(P)  =  f  -  (l/r)  n,- (VO, -VO,)  dS  +  f  (0,  ~02)  n,  V  (l/r)  dS 

K  K 

+  f  A0n,  V  (l/r)  dS  +  4nO„.  (17) 

JSI3 

If  P  is  on  the  surface  of  S12,  the  potential  at  P  can  be  derived  in  a  similar 
manner  by  repeated  application  of  Eq.  8;  this  gives 

2nO,(P)  +  2n02(P)  =  f  -  (l/r)  n,  (VO,  -  V02)  dS  +  f  n,  •  V  (l/r)  (O,  -  02)  dS 

Js,2  Js12 

+  f  AOn,-V(l/r)  dS  +  4ttO_  (18) 

Js13 

where  0((P)  and  02(P)  are  the  potentials  at  point  P  on  the  external  and  internal 
sides  of  Q2  of  the  surface  S12,  respectively.  Equation  18  can  be  rewritten  in  the 
following  two  forms: 

4nO  (P)  =  f  -  (l/r)  n -(VOj-VOj)  dS  +  f  n,  V  (l/r)  (0,-02)  dS 

K  K 

+  f  AOn,-V  (l/r)  dS  +  4nO„  +  2n  [0,(P)-02(P)]t  (19) 

JS .. 
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and 


4i*t>2(P)  = 


(1/r)  ^-(VO.-VOj)  dS  + 


n,  V  (1/r)  (0,-Oj)  dS 


+ 


AOiyV  (1/r)  dS  +  4nO„.  -  2n  [<t>t(P)~ <t>2(P)]. 


(20) 


The  velocity  at  point  P  can  be  obtained  by  applying  the  Vp  operator  (the 
derivatives  being  evaluated  with  respect  to  the  coordinates  of  point  P)  to  Eqs. 
17,  19,  or  20,  depending  on  the  location  of  point  P.  For  illustrative  purposes, 
the  application  of  the  Vp  operator  is  as  follows: 


4nVp<J>1(P)  =  f  Kj  (VO,-V<t>2)  dS  +  f  K4  (<t>f-<t>2)  dS 

K  K 

+  f  K4  AO  dS  +  4nV.  +  2n  [VpO,(P)- Vp02(P)] 

Js., 


=  4nV(P) 


(21) 


where  K,  [=  -Vp(l/r)  =  t/i3]  and  K4  [=  Vp[^,  V(l/r)]  =  ^  V^f/r3]  are 
vector  kernel  functions  which  depend  only  on  the  geometry  and  are  independent 
of  flow  or  boundary  conditions. 

As  a  consequence  of  Green’s  identit;  Eq.  19  and  Eq.  21  show  that  0,(P) 
and  V,(P),  the  potential  and  velocity  of  a  field  point  P  in  the  external  region  Q 
or  on  the  external  surface  of  S,2,  can  be  expressed  in  terms  of  the  distribution 
of  n,  •  (VO  -  V02),  (O,  -  02)  on  S12,  and  AO  on  he  surface  of  the  wake.  The  value 
0,(P)  or  Vj(P)  can  be  uniquely  determined  on  "  if  boundary  conditions  for  02 
on  the  internal  flow  region  Q2  are  specified.  1 .  e  boundary  conditions  for  02 
depend  on  the  physical  problems  encountered.  They  may  also  affect  the 
numerical  error  and  stability  of  the  solution  of  0,(P)  and  V,(P). 

A  point  source  of  strength  a  at  source  point  Q  induces  a  potential  «t>t(P) 
and  velocity  Vf(P)  at  a  point  P7  such  that 


*,( P)  = 


o 

4nr 


and 


V(P)  = 


o  r 
4n  jj 


(22) 


where  r  is  the  magnitude  of  the  vector "r  pointing  from  Q  to  P.  A  doublet  of 
strength  h  at  source  point  Q  induces  a  potential  <t>d(P)  and  velocity  ^d(P)  at 
point  P  such  that 
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and 


(23) 


=  ‘  i;^'V1/r)  *  ~k ^Vi/o 


Vd(P)  =  ipVp[frVp(l/r)] 


where  "n  is  the  normal  along  the  axis  of  the  doublet  pointing  from  the  negative 
to  the  positive  end  of  the  doublet.  The  derivatives  in  the  vector  operator  Vp  and 
VQ  are  taken  with  respect  to  the  coordinates  of  points  P  and  Q,  respectively. 

With  Eq.  22,  the  first  integral  on  the  right-hand  side  of  Eq.  19  can  be 
interpreted  as  the  potential  induced  at  point  P  due  to  a  source  distribution  on 
surface  S,2  whose  strength  is  equal  to  the  boundary  value  of  •  (VO,  -  V<t>2). 

With  Eq.  23,  the  second  integral  of  the  same  equation  can  be  interpreted  as  the 
potential  induced  at  point  P  due  to  a  distribution  of  doublets  on  S,2  whose  axes 
lie  along  with  the  unit  normal  surface  vector  n,  and  whose  strength  is  the 
boundary  value  of  -  (<t>,  -  <t>2).  Note  that  the  derivatives  of  the  vector  operator  A 
in  this  integral  are  taken  with  respect  to  the  coordinates  of  source  point  Q  on 
surface  S,2.  Equation  19  can  now  be  written  as 


4nO,(P)  = 


(K.o-K^i)  dS 


+ 


dS 


+  4n0m  —  2nf/(P)  (24) 

where  K,  [=  -  (1/r)]  and  K,  [=  n,  Vp(l/r)  *  n", -t/V]  are  scalar  kernel 
functions  which  depend  only  on  the  geometry  of  the  boundary  and  wake  shape 
and  are  independent  of  the  flow  and  boundary  conditions,  and  pw  is  the  doublet 
strength  at  the  wake  surface  to  be  determined  from  a  Kutta  condition.  Likewise, 
Eq.  21  can  be  written  as 

4nV(P)  =  f 

J  S 

where  nT,  is  the  outward  normal  at  point  P. 

The  concept  of  placing  singularities  such  as  sources  and  doublets  of 
specified  strength  over  the  boundary  surfaces  of  a  flow  field  forms  the  basis  of 
the  surface  singularity  technique.  As  discussed  previously,  various  assumptions 
and  boundary  conditions  regarding  the  internal  field  may  be  used  to  obtain  a 
unique  solution  for  problems  which  are  physically  meaningful.  Different 
assumptions  demand  different  solution  techniques.  The  choice  of  assumption 
depends  on  the  nature  of  the  problem,  such  as  lifting  or  nonlifting  bodies,  thick 
or  thin  wing  sections,  etc.  If  the  integral  equations  are  to  be  solved  numerically, 
the  assumptions  play  an  important  role  in  the  errors  and  accuracy  of  the 
solutions.  Some  interesting  and  practical  formulations  are  given  in  the  following 
sections. 


K30  dS  -  f  K4p  dS  +  f  Kj*  dS  +  4tiV„  +  2no(P)nj  (25) 

,2  jS,2  Js„ 
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SOURCE-ONLY  FORMULATION 

This  formulation  was  adopted  by  Hess  and  Smith2  to  solve  the  potential 
flow  about  arbitrary  nonlifting  bodies.  The  internal  potential  <t>2  is  assumed  to 
be  continuous,  and  at  the  wetted  surface  it  has  the  same  value  as  the  external 
total  potential  (4* ,  —  <t>2  =  0  at  Si2  in  Fig.  1).  The  normal  derivative  of  0  is 
discontinuous  at  the  surface.  After  the  inner  product  of  Eq.  25  and  the  normal 
vector  nt  at  point  P  has  been  taken  and  the  terms  grouped,  Eq.  25  becomes 

2no(P)  +  f  (Kj-n^dS  =  ~4n  (V.- V(P))  (26) 

Js12 

Equation  26  is  an  integral  equation  since  the  unknown  source  strength  o  appears 
under  an  integral  sign.  The  equation  is  identified  as  the  Fredholm  integral 
equation  of  the  second  kind.  Numerical  discretization  converts  the  integral 
equation  into  a  system  of  linear  equations  denoted  by 

[A]o  =  b 

where  [A]  is  the  influence  coefficient  matrix,  and  b  is  the  boundary  condition,  a 
is  to  be  solved  at  the  discretized  points.  With  the  presence  of  the  singular  term 
2no(P)  in  Eq.  26,  [A]  becomes  diagonal  dominating  and  well  conditioned. 

DOUBLET-ONLY  FORMULATION 

The  normal  gradients  of  <t>,  and  <t>2  are  assumed  to  be  equal  at  the  wetted 
surface,  separating  the  region  Q,  from  the  region  Q2  (n:(V<t>,-V<l>2)  =  0  at  SJ2). 
After  implementation  of  this  assumption,  the  o  terms  vanish  from  Eqs.  24  and 
25.  This  leads  to  the  doublet-only  formulation.  Physically,  this  assumption 
implies  that  the  normal  component  of  the  velocity  at  the  inner  and  outer 
surfaces  of  the  wetted  surface  S,2  is  continuous.  If  the  wetted  surface  is  solid, 
then  Vn=0  at  the  inner  surface  and  <t>2  becomes  constant  everywhere  in  Q2. 

After  the  o  term  has  been  dropped  and  the  inner  product  with  the  normal  vector 
li,  taken,  Eq.  25  becomes 

4nV(P)-'nJ  =  4nV.-n,  -  f  M(K4  n,)  dS  +  f  Mw(K4  n,)  dS.  (27) 

Jsl2  K 

Equation  27  is  a  Fredholm  integral  equation  of  the  first  kind  with  p  to  be 
solved.  Through  numerical  discretization,  Eq.  27  can  be  converted  to  a  system 
of  linear  equations  [A]^  =  b  where  [A]  is  singular;  that  is,  the  sum  of  the 
elements'  in  every  row  is  zero.  The  system  of  equations  is  not  linearly 
independent.  To  obtain  a  unique  solution,  an  additional  condition  is  required. 
One  of  the  conditions  is  to  specify  a  fixed  value  of  p  at  a  given  point  on  S|2. 

SOURCE  AND  DOUBLET  FORMULATION 

Equation  20  expresses  the  potential  at  point  P  on  surface  S,2  on  the  Q2  side 
and  can  be  rewritten  in  terms  of  source  and  doublet  distributions  as 
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4n<t>2(P)  = 


(K.o.-K^)  dS 


dS  +  2n^i(P)  +  4710,,, 


(28) 


where  K,  and  Kj  are  the  scalar  kernel  functions  as  defined  in  Eq.  24.  The 
obvious  external  boundary  condition  on  the  wetted  surface  is 


n.-VO,  =  Vn.  (29) 

Vn  is  zero  if  the  surface  is  solid.  As  mentioned  earlier,  to  obtain  a  unique 
solution  of  Eq.  28,  some  assumption  regarding  the  internal  flow  field  should  be 
made.  Johnson,®  Bristow,9  and  Maskew5, 10  have  applied  a  Dirichlet  condition  to 
the  internal  flow  field  and  set  0,  =  <t>  .  With  these  defined  internal  and  external 
boundary  conditions,  the  source  strength  on  the  surface  can  be  written 

c  =  n,  •  (VO,  -  V02)  =  Vn-n,  •  V..  (30) 

Suppose  the  lifting  body  under  consideration  is  a  propeller  blade  which 
rotates  about  an  axis  with  angular  velocity  co.  The  analysis  derived  previously 
still  holds  with  respect  to  a  moving  frame  fixed  to  the  propeller  blade.  In  this 
moving  frame  fixed  to  the  propeller,  the  boundary  condition  expressed  in  Eq.  29 
becomes 

=  V-Dj'SxR  (31) 

where  R  is  the  position  vector  of  point  P  on  the  surface  of  the  propeller  with 
respect  to  the  frame  fixed  to  the  propeller.  Equation  30  becomes 

o  =  n,  •  (VO,  -  V4>2)  =  Vn-n,  -V^-n,  -wx  R.  (32) 

All  quantities  on  the  right-hand  side  of  Eq.  32  are  known.  It  is  now  clear  that 
application  of  the  Dirichlet  boundary  condition  on  the  internal  region 
determines  the  source  strength  o  in  Eq.  28.  Equation  28  can  be  written  as 

2nM(P)  = 

and  Q  is  a  source  point  on  surface  S|2  excluding  P.  Equation  28  is  a  Fredholm 
integral  equation  of  the  second  kind  and  can  be  solved  uniquely  in  terms  of  the 
doublet  strength  p. 


K.o  dS 
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-  f  KjM.  dS  +  f  K^Q)  dS  (33) 


SOURCE  AND  VORTICITY  FORMULATION 

A  constant-strength  doublet  panel  is  equivalent  to  a  ring  of  line  vortices 
around  the  perimeter.  Hess11  proved  in  a  more  general  case  that  a  surface 
doublet  distribution  of  density  p  can  be  replaced  by  a  vortex  sheet  (on  the  same 
surface  as  the  doublet  sheet)  plus  a  concentrated  vortex  filament  around  the 
edge  of  the  sheet.  The  vorticity  "y  on  the  sheet  can  be  related  to  the  doublet  as 
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y  =  nxVfi 


where  n  is  the  unit  normal  vector  of  the  surface  and  Vp  is  the  gradient  of  the 
doublet  strength.  The  strength  of  the  concentrated  vortex  filament  around  the 
edge  is  equal  to  the  local  edge  doublet  strength. 

In  the  three-dimensional  case,  the  magnitude  of  the  strength  and  the 
direction  of  the  vortex  are  unknown.  Care  must  be  taken  to  ensure  that  Kelvin’s 
circulation  condition  (zero  divergence)  is  satisfied  everywhere.  The  numerical 
implementation  of  the  vorticity  formulation  is  more  complicated  than  for  the 
source/doublet  formulation.  Bristow9  preferred  the  source/doublet  formulation 
over  the  source/vorticity  formulation  for  the  design  problem. 

Although  all  types  of  singularity  distributions  just  discussed  can  be  derived 
directly  or  indirectly  from  Green’s  identity,  the  resulting  numerical  formulations 
and  accuracies  are  quite  different.  Bristow9  showed  that  the  source/doubiet 
distribution  in  general  provides  a  source  distribution  milder  than  the  source-only 
solution  and  a  doublet  distribution  milder  than  the  doublet-only  distribution.  In 
summary,  the  source/doublet  distribution  is  attractive  because 

1 .  Relatively  mild  singularity  distributions  suppress  the  numerical 
instabilities,  which  may  otherwise  be  prevalent  for  these  high  lift  geometries. 

2.  Direct  relationships  between  velocity  and  singularity  strengths  on  the 
boundary  surface  simplify  calculations;  which  facilitates  implementation  of  the 
boundary  layer  displacement  formulation  for  inviscid/viscous  interaction. 

3.  Singularity  distributions  vanish  as  the  perturbation  field  vanishes, 
thereby  eliminating  possible  residual  errors. 

4.  With  a  source/doublet  distribution,  it  is  easier  to  obtain  a  numerically 
stable  solution  in  the  design  problem. 

NUMERICAL  DISCRETIZATION 

The  conversion  of  a  differential  equation  to  an  integral  equation,  illustrated 
in  a  previous  section,  becomes  a  major  technique  for  solving  initial-value  and 
boundary-value  problems  of  ordinary  and  partial  differential  equations.  Only  in 
a  limited  number  of  cases  can  the  Fredholm  equation  be  solved  in  closed 
analytical  forms.  In  general,  these  equations  must  be  solved  numerically.  The 
Fredholm  integral  equation  of  the  first  kind  is  more  difficult  to  solve  than  the 
second  kind.  For  this  reason,  the  discretization  of  Eq.  33  with  appropriate 
boundary  conditions  will  be  addressed  here.  The  complete  discretization  of  the 
problem  involves  three  different  tasks;  (a)  discretization  of  the  boundary 
condition,  (b)  discretization  of  surface  geometry,  and  (c)  discretization  of  the 
singularity  distribution. 

DISCRETIZATION  OF  THE  BOUNDARY  CONDITION 

A  continuous  solution  of  an  integral  equation  demands  that  the  boundary 
condition  be  satisfied  on  the  physical  boundary  in  a  continuous  manner.  With  a 
numerical  approach,  the  boundary  condition  can  be  satisfied  only  at  a  finite 
number  of  selective  collocation  points.  Consequently,  the  velocity  components 
and  potential  at  the  physical  surface  between  the  collocation  points  are  not 
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likely  to  satisfy  the  imposed  condition.  For  example,  if  the  normal  component 
Vn  of  the  velocity  is  prescribed  at  the  boundary,  the  discretized  solution  will  give 
the  potential  field  that  produces  the  correct  velocity  component  only  at  the  discretized 
collocation  points  on  the  boundary.  Between  the  collocation  points,  the  calculated 
potential  field  in  general  will  give  a  value  for  VD  that  is  quite  different  from  the 
physical  one.  This  difference  is  sometimes  referred  to  as  leakage.  The  problem  is 
inherent  in  discretization.  The  leakage  may  not  be  eliminated  completely  but  can 
be  reduced  by  using  a  more  accurate  surface  discretization  (more  collocation 
points  when  the  surface  normal  vector  varies  rapidly  and  a  smoother  distribution 
of  singularity,  such  as  provided  by  a  higher  order  model). 

DISCRETIZATION  OF  THE  GEOMETRICAL  SURFACE 

The  smooth  continuous  shape  of  the  body  surface  is  represented  by  a  number 
of  plane  quadrilateral  panels  whose  corners  are  the  projection  of  the  surface  points 
on  the  panels.  Use  of  a  flat  panel  to  approximate  a  curved  surface,  makes  it 
inevitable  that  geometrical  discontinuities  will  exist  on  the  surface  of  the  panels 
where  the  singularities  are  distributed.  To  approximate  the  continuous  physical 
surface  meaningfully  with  a  collection  of  flat  panels,  the  deviations  of  the 
surface  points  and  their  projection  on  the  panel  plane  should  be  kept  as  small  as 
possible.  On  a  surface  area  where  curvature  is  large,  small  panels  should  be 
used.  The  centroid  of  the  panel  plane  can  be  taken  as  the  collocation  point  on 
which  the  boundary  conditions  are  applied. 

The  collocation  point  on  the  centroid  of  the  plane  panel  may  not  be  on  the  physical 
surface,  and  the  unit  normal  vector  of  the  panel  may  not  be  aligned  with  the  normal 
of  the  physical  surface.  The  deviation  between  the  true  and  numerical  solution  depends 
on  the  distribution  and  orientations  of  the  panels.  It  should  be  kept  in  mind 
that  the  numerical  solution  is  meaningful  only  at  the  collocation  points. 

DISCRETIZATION  OF  SINGULARITY  DISTRIBUTION 

The  singularities  described  previously  are  distributed  on  the  panel  surface 
rather  than  on  the  physical  surface.  The  strengths  of  the  singularities  can  be 
constant  or  can  vary  continuously  over  the  panel,  depending  on  whether  first-  or 
higher  order  approximations  are  used.  For  a  simple  singularity  distribution,  the 
influence  of  a  single  panel  of  arbitrary  shape  at  a  given  field  point  can  be 
computed  entirely  analytically.  Some  examples  can  be  found  in  Hess  and  Smith1 
and  Newman.*  Hie  discontinuity  of  the  source  strength  at  the  edge  a  of  panel 
may  result  in  an  erroneous  velocity  which  becomes  logarithmically  unbounded  in 
the  vicinity  of  the  edge,  and  the  discontinuity  of  the  doublet  strength  may  cause 
a  jump  of  potential  and  velocity  near  the  edge.  To  reduce  the  error,  a  consistent  higher 
order  approximation  is  required  to  provide  continuous  source  and  doublet  strengths 
on  some  points  along  the  edge.  The  discontinuities  of  singularity  strength  may  also 
be  eliminated  by  weighted  averaging  along  a  common  panel  edge.  The  improved 
accuracy  is  obtained  at  the  expense  of  computational  effort.  The  accuracy  of  a 
first-order  approximation  may  be  improved  by  using  more  panels  on  surface 
regions  where  curvature  changes  rapidly.  The  choice  of  a  first-  or  higher  order 
approximation  should  be  decided  by  balancing  cost  and  accuracy. 

•Newman,  J.N.,  “Distribution  of  Sources  and  Normal  Dipoles  Over  a  Quadrilateral  Panel,” 
(submitted  for  publication,  1986). 
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NUMERICAL  SOLUTION 

Numerical  discretization  and  solution  of  the  integral  equation,  Eq.  33,  are 
discussed  here.  For  illustration  purposes,  the  singularities  (source/doublet)  are 
assumed  to  have  constant  strength  on  a  given  panel.  Suppose  that  the  physical 
surfaces  are  replaced  with  N  flat  panels  and  the  wake  surface  is  replaced  with  M 
panels.  Equation  33  can  then  be  approximated  as 

N  N  M 

2W.  =  ^  K^,  dS  -  K,o.  dS  -  ^  OO,  ^  ^  <34> 

i=  1  S.  i=  1  S,  i=  1  S. 

l*j 

where  p.  is  the  doublet  strength  on  the  panel,  S.  is  the  area  of  the  panel,  and  o.  is 
the  source  strength.  In  the  current  example,  a  is  known  and  assumed  to  be  constant 
over  the  panel  area  S..  Since  the  strengths  of  the  singularities  are  constant  on  each 
panel,  they  can  be  taken  out  of  the  integral  signs  in  Eq.  34.  The  scale  kernel  K, 
and  Kj  can  be  integrated  over  the  panel  surface  analytically.1'*  Equation  34  can 
then  be  written  in  a  simplified  form  as 

N  N 

Inn.  =  *  C«  ~  ]E  °<  B« 

i=  1  i=  1 

l*j 

where  B  and  C|]  matrices  are  the  surface  integration  of  the  kernel  K,  and  over 
the  panel  containing  field  point  i  with  respect  to  the  field  point  j.  Physically, 
they  express  the  potential  induced  at  collocation  point  j  by  the  source  and  doublet 
distribution  on  the  panel,  which  contains  collocation  point  i.  When  the  collocation 
point  j  is  sufficiently  far  away  from  point  i,  the  full  expression  is  replaced  by  a 
simpler  formulation  which  yields  results  within  a  reasonable  tolerance.1  The 
doublet  strength  on  the  wake  surface  is  determined  from  a  Kutta  condition. 
Equation  35  can  be  written  in  matrix  form 

[A]  X  =  B  (36) 

where  [A]  is  the  influence  coefficient  matrix  which  has  a  dimension  of  N  by  N,  X 
is  the  doublet  strength  column  vector  with  dimension  N,  and  B  is  the  boundary 
condition  vector  with  dimension  N.  The  singular  term  of  the  Fredholm  integral 
equation  of  the  second  kind  makes  the  [A]  matrix  well  conditioned  (in  general, 
it  is  diagonal  dominating). 

In  actual  numerical  computation,  the  assembly  of  the  matrix  [A]  involves 
the  evaluation  of  By  and  Ctj  and  is  very  time  consuming.  [A]  is  a  full  matrix, 
and  the  solution  of  Eq.  36  demands  also  a  large  amount  of  computing  effort, 
especially  when  the  number  of  panels  is  large.  The  block  Gauss-Seidel  method 
can  be  used  in  solving  such  a  large  system. 

•Newman,  J.N.,  “Distribution  of  Sources  and  Normal  Dipoles  Over  a  Quadrilateral  Panel,” 
(submitted  for  publication,  1986). 
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VISCOUS/INVISCID  INTERACTION 

In  virtually  every  aspect  of  hydrodynamic  problems  the  viscous  effects  are 
significant  and  must  be  accounted  for.  Although  the  Navier-Stokes  equation 
applies  to  such  flows,  it  requires  significant  computational  resources7,  especially 
at  high  Reynolds  number  conditions.  For  this  reason,  simpler  viscous/inviscid 
interactive  approaches  were  derived. 

In  a  recent  review  paper  on  calculation  methods  for  flow  on  airfoils, 

Cebeci,  Stewartson,  and  Whitelaw12  concluded  that  the  interactive  approach  is 
more  economical  than  the  direct  application  of  Reynolds-averaged  two- 
dimensional  forms  of  the  Navier-Stokes  equations.  In  the  interactive  approach, 
the  flow  is  divided  into  an  outer  inviscid  region  and  a  thin  viscous  region  near 
the  solid  surface  and  around  the  wakes.  The  inviscid  region  can  be 
approximated  by  panel  methods  if  the  flow  is  irrotational.  The  interaction 
between  inner  and  outer  regions  can  be  established  by  two  different  methods. 
The  first  method  as  presented  by  Mahgoub  and  Bradshaw,13  is  differential.  The 
regions  are  matched  at  an  arbitrary  line  that  encloses  the  solid  surface.  The 
inner  solution  is  obtained  by  finite-differencing  the  parabolized  Navier-Stokes 
equation,  and  the  outer  solution  is  obtained  by  the  panel  method.  The  second 
method  is  integral.  In  this  method,  the  inner  region  is  not  physically  prescribed, 
and  the  outer  region  contacts  the  solid  surface  directly.  The  viscous  effect  is 
extended  to  the  inviscid  region  by  the  deficit  formulation,  as  described  by 
LeBalleur.14  Both  differential  and  integral  methods  can  able  to  calculate  the  flow 
patterns  beyond  the  separation  points  with  reasonable  accuracy.  Both  methods 
can  be  extended  to  lifting  body  calculations  in  which  the  wake  curvature  effects 
are  also  important.  Although  the  differential  method  gives  more  detailed 
information  about  the  flow  field,  it  is  impractical  for  design  purposes  because  of 
the  huge  computational  resource  required.  Kline,  Cantwell,  and  LilleyiS  showed 
that  an  integral  boundary  layer  method  need  not  be  less  accurate  than  methods 
based  on  Reynolds  averaged  equations. 

There  are  two  popular  models  of  the  viscous  displacement  effect  in  integral 
methods.  One  model  requires  the  computation  of  boundary  layer  displacement 
thickness,  and  the  zero-normal-velocity  boundary  condition  is  then  applied  to 
the  modified  geometry  as  described  by  Dutt  and  Sreekanth.16  The  other  model 
was  originally  proposed  by  Lighthill.17  It  applies  a  relatively  simple  transpiration 
condition  on  the  physical  surface.  The  latter  model  is  more  practical  if  a  panel 
method  is  being  used  to  compute  the  inviscid  flow,  since  the  matrix  of  influence 
coefficients  of  Eq.  36  need  not  be  computed  during  the  interactive  process  once 
it  is  set  up.  On  the  basis  of  the  boundary  transpiration  concept.  Green1*  derived 
an  entrainment  method  to  predict  the  flow  in  a  turbulent  boundary  layer.  This 
method  has  been  further  improved  by  East,19  Lock  and  Firnin,20  Le  Balleur,14 
and  Melnik  and  Brook21  to  account  for  flow  separation.  Other  approaches 
capable  of  treating  flow  separation  were  given  by  Dvorak,  Maskew,  and 
Woodward,22  Dvorak,  Woodward,  and  Maskew,23  and  Maskew,  Rao,  and 
Dvorak.24  Application  of  a  panel  method  to  attached  cavitation  has  been 
reported  by  Franc  and  Michel.25 
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APPLICATION  OF  PANEL  METHOD  AS  A  DESIGN  TOOL 

One  of  the  design-type  problems  involves  the  determination  of  body 
contour  that  will  produce  a  prescribed  pressure  distribution.  The  panel  methods 
have  been  successful  in  predicting  flow  pattern  and  surface  pressure  distribution 
involving  complicated  geometries  in  potential  flow,  but  they  are  not  widely  used 
in  designing  lifting-surface  elements.  One  of  the  reasons  is  that  the  design 
problem  is  more  complicated  than  the  analysis  problem,  and  the  design 
formulation  is  often  unstable  unless  some  precautions  are  taken.  The  other 
reason  is  that,  due  to  the  non-linearity  of  the  boundary  conditions,  an  iterative 
approach  is  needed  to  obtain  a  solution  and  the  computing  cost  becomes 
prohibitive.  Several  design  methods  based  on  the  panel  approach  have  been 
reported  in  the  past  for  designing  wing  section  contours  that  will  produce  a 
prescribed  pressure  distribution,26,27  but  not  all  are  satisfactory.  Some  of  the 
drawbacks  are  (a)  the  calculations  fail  to  coverage,  (b)  the  design  contour  is 
unrealistically  wavy,  and  (c)  computing  cost  is  prohibitive.  Slooff28  has  discussed 
several  other  design  methods. 

Bristow29  and  Hawk30  developed  a  perturbation  analysis/design  method  for 
analyzing  a  series  of  arbitrary  small-geometry  perturbations  to  a  baseline 
configuration  and  for  designing  a  wing  section  that  will  produce  a  prescribed 
pressure  distribution.  The  method  is  based  on  the  panel  method  with 
source/doublet  surface  singularity  distributions  described  previously  (constant 
source  and  quadratic  doublet  distribution).  The  authors  reported  that  the 
method  is  cost  effective,  accurate,  and  stable. 

The  perturbation  analysis  method  was  developed  for  efficient  and  accurate 
analysis  of  a  series  of  arbitrary  small  geometry  perturbations  to  a  baseline 
configuration.  The  velocity  potential  for  the  perturbed  configuration  is  obtained 
by  a  first-order  Taylor  series  expansion  about  the  velocity  potential  of  the 
baseline  configuration.  Basically,  it  is  a  linear  extrapolation  procedure  and 
bypasses  the  two  most  expensive  steps:  (a)  assembling  the  new  influence 
coefficient  matrix,  and  (b)  solving  a  large  system  of  linear  algebraic  equations 
for  the  velocity  potential.  The  authors  attributed  the  success  of  the  method  to 
the  fact  that  only  the  velocity  potential — not  velocity  and  pressure — is  linearized 
with  respect  to  the  geometry  perturbation  coordinate,  and  the  nonlinear  terms 
are  much  smaller  for  velocity  potential  than  for  either  velocity  or  pressure.  The 
method  is  applicable  to  large  perturbations  of  wing  thickness,  camber,  and  twist. 

The  perturbation  design  method  is  the  logical  extension  of  the  perturbation 
analysis  method  already  described.  The  design  can  be  initiated  from  a  baseline 
geometry  obtained  from  lifting  surface  theory.  The  geometry  is  perturbed  in 
such  a  manner  that  the  pressure  distribution  on  the  perturbed  surface  satisfies  a 
“target”  distribution.  Iteration  is  required  because  the  pressure  is  a  nonlinear 
function  of  the  geometry  perturbation.  In  each  iteration  cycle,  there  are  two 
major  steps:  (a)  an  analysis  solution  of  the  geometry;  in  the  first  iteration,  the 
analysis  is  performed  on  the  baseline  configuration;  and  (b)  a  modification  of 
the  geometry.  The  first  step  follows  the  analysis  procedure  already  described. 

The  modification  of  geometry  in  the  second  step  is  obtained,  for  example,  by 
requiring  that  the  pressure  distribution  on  the  modified  contour  satisfies  a 
“target”  distribution.  The  first-order  Taylor  series  expansion  of  the  velocity 
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potential  about  the  baseline  geometry  requires  knowledge  of  the  partial 
derivatives  of  velocity  potential  with  respect  to  the  geometry  coordinates.  The 
partial  derivatives  are  calculated  once  only  and  then  are  stored  for  repetitive  use  during 
the  iteration  procedure.  Hawk  and  Bristow30  reported  that  the  perturbation  analysis 
method  is  competitive  in  accuracy  with  that  of  conventional  panel  analysis 
methods,  the  computing  cost  of  each  successive  application  is  one  order  of 
magnitude  less,  and  the  perturbation  design  method  is  efficient  in  calculating 
three-dimensional  wing  section  geometry  corresponding  to  a  prescribed  pressure 
distribution.  Although  the  sample  problems  presented  by  Bristow29  and  Hawk30 
are  wing-section  related,  the  extension  of  the  method  to  propeller  blade  design  is 
straightforward. 


SAMPLE  CALCULATIONS 

Lifting  surface  theory  and  equivalent  two-dimensional  theory  have  been 
used31,32  to  predict  the  pressure  distribution  on  the  surface  of  propeller  blades. 
Neither  theory  accounts  for  the  effect  of  the  hub.  Due  to  this  limitation  in 
design  and  analysis  procedures,  the  hub  effect  is  often  greatly  simplified  or 
ignored.  Detailed  knowledge  of  the  pressure  distribution  on  the  blade  surface— 
especially  around  the  leading  edge — is  important  to  cavitation  inception 
prediction.  However,  neither  theory  performed  satisfactorily,  especially  at  the 
off-design  conditions. 

The  panel  method  distributes  the  singularities  on  the  surfaces  of  both  the 
blade  and  the  hub.  A  hub  with  any  arbitrary  shape  can  be  easily  modeled  as  an 
integrated  part  of  a  propeller.  The  geometry  of  the  propeller  unit  can  be 
modeled  in  as  much  detail  as  the  maximum  number  of  panels  allows,  avoiding 
the  shortcoming  suffered  by  the  lifting  surface  theory. 

Benchmark  comparisons  for  the  panel  method  calculations  presented  in  this 
report  were  performed  using  the  experimental  results  previously  reported  at  the 
David  Taylor  Research  Center  (DTRC).  Model  tests  were  conducted  with  DTRC 
Propeller  4718  with  a  hub/diameter  ratio  of  0.3.  The  propeller  drawing  is  shown 
in  Fig.  3,  with  tabulated  geometry  listed  in  Table  1.  The  propeller  was  tested  in 
an  open  jet  test  section  of  the  DTRC  36-in.  variable-pressure  water  tunnel  and 
in  open  water  in  the  DTRC  high-speed  towing  basin.  The  measurements  of  the 
pressure  on  the  blade  surface  were  taken  at  locations  that  were  0.5,  0.7,  0.8, 
and  0.9  of  the  tip  radius. 

During  the  experiment,33  an  extensive  calibration  program  was  conducted  to 
arrive  at  accurate  calibrations  for  test  purposes,  and  to  investigate  possible 
systematic  errors  in  pressure  measurement  instrumentation.  For  any  single 
calibration,  the  error  band,  based  on  a  95%  confidence  level,  was  calculated 
from  the  standard  deviation  relative  to  a  straight  line  calculated  sensitivity.  It 
was  determined  that  the  uncertainty  of  pressure  measurements  at  low  test  speed 
(6  knots)  and  at  high  test  speed  (1 1  knots)  is  ±  0.05  psi  and  ±  0.07  psi, 
respectively.  The  detailed  method  of  calibration  can  be  found  in  Reference  33. 

The  panel  method34"36  used  in  the  calculations  was  the  VSAERO  code;  it  has 
piecewise  constant  source  and  doublet  formation  as  described  previously.  The  exact 
blade  surface  was  approximated  with  609  flat  quadrilateral  panels,  and  the  hub 
was  approximated  with  144  panels.  The  panelized  propeller  is  shown  in  Fig.  4. 
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To  demonstrate  the  hub  effect,  the  performances  of  the  propeller  without  and 
with  a  hub  were  calculated.  The  results  were  then  compared  with  predictions 
from  lifting  surface  theory  and  with  experimental  data. 

Comparisons  were  made  between  the  measured  surface  pressure  and 
predictions  from  the  lifting  surface  theory  and  from  the  panel  method  without 
the  hub.  The  results  are  shown  in  Fig.  5.  In  general,  predictions  using  the  panel 
method  were  an  improvement  over  the  lifting  surface  theory.  However, 
predictions  at  the  0.5  radius  were  less  than  satisfactory.  The  large  discrepancies 
between  the  predictions  and  experimental  data  were  attributed  to  the  hub  effect 
that  was  not  accounted  for  in  the  calculations. 

With  the  hub  model  included,  the  pressure  distributions  were  again 
calculated  by  the  panel  method,  and  the  results  are  shown  in  Fig.  6.  The 
agreement  between  the  predictions  and  experimental  data  improved  significantly, 
especially  at  the  0.5  radial  section  where  the  hub  effect  is  expected  to  be  great. 

The  fluid  accelerates  as  it  passes  the  hub  and  influences  the  flow  at  the 
blade  sections.  Locally,  each  blade  section  operates  at  an  advance  ratio  that  is 
different  from,  and  greater  than,  the  design  value.  The  difference  is  greatest 
near  the  hub/blade  junction  and  diminishes  as  the  ratio  r/R  increases. 
Qualitatively,  the  hub  increases  the  pressure  near  the  leading  edge  and  decreases 
the  pressure  near  the  trailing  edge  on  the  suction  side  and  imposes  the  opposite 
effect  on  the  pressure  side.  The  panel  method  predicts  this  effect  correctly,  and 
the  agreement  with  the  experiment  data  is  excellent. 

Predictions  of  pressure  distributions  at  off-design  conditions  were 
calculated,  and  comparisons  with  predictions  from  lifting  surface  theory  and 
experimental  data  are  shown  in  Figs.  7-9.  The  calculated  pressure  peaks  at  the 
leading  edge  agree  well  with  the  experimental  data. 

CONCLUSIONS 

It  is  demonstrated  here  that  the  panel  method  is  an  improvement  over 
lifting  surface  theory  in  predicting  pressure  distribution  on  propeller  blade  surface. 
The  most  significant  improvements  are  that  (a)  the  panel  method  can  accommodate 
complicated  hub  geometry  with  relative  ease  and  (b)  it  can  provide  more  reliable 
information  near  the  leading  edge  where  many  other  approaches  have  failed, 
especially  under  off-design  conditions.  Although  the  formation  of  the  panel  code 
VSAERO  is  first  order  in  nature,  the  results  agree  with  the  experimental  data 
well.  One  computation  (for  a  given  advance  coefficient  J)  requires  23  minutes 
CPU  on  a  VAX/780  machine. 
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FRACTION  OF  CHORD  FROM  LEADING  EDGE,  xc 
Fig.  5b.  xR  =  0.7. 

Fig.  5.  Comparison  of  blade  pressure  measurement  Cp  with  lifting  surface 
theory  and  panel  method  calculations  (without  hub  model)  for 
Propeller  4718  at  various  fractions  of  propeller  radius  x„. 


Fig.  5c.  xR  =  0.8. 


Fig.  5d.  xR  =  0.9. 
Fig.  5.  (Continued) 
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Fig.  6a.  x„  =  0.5. 
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Fig.  6b.  xR  =  0.7. 

Fig  6.  Comparison  of  blade  pressure  measurements  C„  with  lifting  surface 
theory  and  panel  method  calculations  (with  hub  model)  for 
Propeller  4718  at  various  fractions  of  propeller  radius  xR. 
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Fig.  8d.  xc  =  0.7. 
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Table  1.  Properties  of  Propeller  4718. 

Blade  Geometry  Et  Design 


Diameter  D  =  2.000  ft  (0.610  m) 

Blade  Thickness  Fraction  =  0.069 

Rotation  =  Right  Hand 

Design  Advance  Coefficient  J  =  0.751 

Number  of  Blades  Z  =  3 

Design  Thrust  Loading  Coefficient  CTh  = 

=  0.248 

Hub/Diameter  Ratio  Dh/D  = 

0.30 

Design  Thrust  Coefficient  KT 

=  0.055 

Expanded  Area  Ratio  =0.44 

Design  Torque  Coefficient  KQ 

=  0.0106 

Blade  Parameters 

r/R 

c/D 

P/D 

8t  (deg) 

iQ/D 

t/c 

t/D 

Wc 

VD 

EH 

0.187 

0.718 

-1.65 

0.0 

0.2497 

0.0467 

0.0 

0.0 

0.249 

0.796 

-4.05 

0.0 

0.1771 

0.0441 

0.0044 

0.0011 

0.311 

0.855 

-5.00 

0.0 

0.1280 

0.0398 

0.0085 

0.0027 

19 

0.366 

0.886 

-3.50 

0.0 

0.0910 

0.0333 

0.0099 

0.0036 

0.403 

977^9 

0.40 

0.0 

0.0254 

0.0101 

0.0041 

0.8 

0.409 

0.870 

5.75 

0.0 

0.0469 

0.0192 

0.0097 

0.0090 

MiK:  Mil 

0.365 

0.825 

12.40 

0.0 

0.0419 

0.0153 

0.0082 

0.0030 

0.95 

0.311 

0.786 

16.10 

0.0 

0.0418 

0.0130 

0.0065 

0.0020 

1.0 

0.070 

0.734 

20.00 

0.0 

0.0414 

0.0029 

0.0090 

0.0006 

Thickness  Et  Camber  Distribution 

*c 

!  ET/f 

W 

0.0000 

0.0000 

0.0000 

0.005 

0.0665 

0.0423 

0.0075 

0.0812 

0.0595 

0.0125 

0.1044 

0.0907 

0.025 

0.1466 

0.1586 

0.05 

0.2712 

0.075 

0.2525 

0.3657 

0.1 

0.2907 

0.4482 

0.15 

0.3521 

0.5869 

0.4000 

0.25 

0.4363 

0.7905 

0.3 

0.4367 

0.8635 

0.35 

0.4832 

0.9202 

0.4 

0.4952 

0.9615 

0.45 

0.5 

0.9881 

1.0 

0.55 

0.4846 

0.9971 

0.6 

0.4653 

0.9786 

0.65 

0.4383 

0.9434 

0.4035 

0.75 

0.3612 

0.8121 

0.8 

0.3110 

0.7027 

0  85 

0.2532 

0.5425 

0.1877 

0.3586 

0.95 

0.1143 

0.1713 

0.975 

0.0748 

1.0 

0.0333 

•  NACA  66  section  (DTNSRDC  modified). 

|  t  NACA  a  =  0.8  meanline:  the  design  procedure  determines  the  magnitude  of  the 

camber  at  each  radius  and  uses  the  two-dimensional  chordwise  distribution  of  camber. 
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APPENDIX 

A  fluid  is  incompressible  if  its  particles  maintain  their  density  along  their 
paths,  i.e.,  the  substantial  derivative  of  mass  density  q  is  zero: 


(A.l) 


The  principle  of  mass  conservation  requires  that  the  net  amount  of  mass  flow 
into  a  control  volume  per  unit  time  be  equal  to  the  rate  at  which  the  mass  in  the 
control  volume  is  increasing.  Thus 

|f  +  Vqu  =  0.  (A. 2) 

Equation  A.l  is  the  differential  equation  of  continuity.  The  bold  type  denotes  a  vector 
quantity  J.  From  Eqs.  A.1  and  A.2  it  follows  that  for  incompressible  fluids  the 
equation  of  continuity  is  simply 


Vu  =  0, 


(A.3) 


Whether  or  not  the  flow  is  steady  and  whether  or  not  the  fluid  is  homogeneous. 
Furthermore,  if  the  flow  is  irrotational,  the  circulation  around  a  closed  circuit  is 
zero, 

Ju-dx  =  0.  (A. 4) 

Therefore,  u  dx  is  an  exact  differential,  which  can  be  denoted  by  d<J>,  thus 

a  =  V4>.  (A.5) 

Equations  A.3  and  A.S  imply  that  the  function  <t>  satisfies  the  Laplace  equation, 

V2  <D  =  0.  (A.6) 


Equation  A.6  is  a  kinematic  condition;  velocity  components  can  be  obtained 
from  its  solution.  The  associated  pressure,  however,  can  be  obtained  only  from 
a  dynamic  condition,  that  is,  the  equation  of  motion.  For  propeller  application, 
it  is  more  convenient  tc  express  the  equation  of  motion  with  respect  to  a 
rotating  frame  that  is  fixed  to  the  propeller  axis.  If  the  fluid  is  inviscid  and  the 
reference  frame  is  rotating  with  a  constant  angular  velocity  w  about  the  x  axis, 
Newton’s  second  law  governing  the  flow  becomes 

g  +  Ft,  (a.7) 

where  F  is  the  Coriolis  acceleration  vector,  2  is  the  body  force  potential,  and 
r2  =  y2+z2. 
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The  Coriolis  acceleration  vector  Fc  (0,  -  2cow,  2a>v)  is  perpendicular  to  the 
velocity  vector  (u,  v,  w).  Hence  its  projection  onto  a  streamline  is  zero.  If  the 
flow  is  steady,  then  along  a  streamline  s  the  equation  of  motion  becomes 


If  q  is  constant,  integration  of  Eq.  A. 8  yields 
P  U 2  CO2!3 

—  +  -4-  +  Q - —  =  constant  along  a  streamline.  (A. 9) 

q  2  2 

The  constant  on  the  right-hand  side  of  Eq.  A.9  can  be  determined  by  the 
upstream  condition. 

In  summary,  to  derive  at  Eqs.  A.6  and  A.9,  it  was  assumed  that  (1)  the 
fluid  is  incompressible  and  inviscid,  and  (2)  the  flow  is  irrotational  and  steady. 
As  a  consequence,  the  flow  solutions  can  be  obtained  from  Eqs.  A.6  and  A.9 
instead  of  from  Eqs.  A.2  and  A.7.  (Equations  A.6  and  A.9  are  simplified  forms 
of  Eqs.  A.2  and  A.7.)  Equation  A.6  is  linear  with  linear  boundary  conditions 
(without  free  surface);  it  can  be  solved  very  easily.  The  velocity  components  can 
be  determined  from  Eq.  A.S,  and  the  associated  pressure  calculated  from  Eq. 
A.9.  The  nonlinearity  of  Eq.  A.7  is  reflected  only  in  the  nonlinearity  of  Eq. 

A.9,  and  there  it  presents  no  difficulty  at  all  because  the  nonlinear  term  is 
clearly  determined  and  only  the  pressure  is  to  be  evaluated. 
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